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On the probability distribution of the stochastic saturation scale in QCD 
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It was recently noticed that high-energy scattering processes in QCD have a stochastic nature. An 
event-by-event scattering amplitude is characterised by a saturation scale which is a random variable. 
The statistical ensemble of saturation scales formed with all the events is distributed according to 
a probability law whose cumulants have been recently computed. In this work, we obtain the 
probability distribution from the cumulants. We prove that it can be considered as Gaussian over 
a large domain that we specify and our results are confirmed by numerical simulations. 
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I. INTRODUCTION 



The study of the high-energy limit of QCD has recently received important contributions coming from analogies 
• with reaction-diffusion systems in statistical physics. It has first been shown 1] that the Balitsky-Kovchegov (BK) 
J> I saturation equation Q lies in the same universality class as the Fisher-Kolmogorov-Petrovsky-Piscounov (F-KPP) 
Cn I equation The BK equation describes the evolution with rapidity Y — log(s) of the dipole scattering ampHtude 

■ T(r, Y) where r the dipole size. From the analogy with the F-KPP equation, one can infer that asymptotic solutions 
: of the BK equation are travelling waves. This property states that T{r, Y) is a function of the single variable rQs{Y) 
' where Qs{Y) is the saturation scale. 

\^ , It has then been realized 0,0 that one also has to take into account gluon-numher fluctuations. The resulting set 

■ of equations can be seen ^ as a reaction-diffusion problem. Alternatively, one can write the evolution equation as 
a Langevin equation which lies in the same universality class as the stochastic F-KPP (sF-KPP) equation 5]. This 
is a Langevin equation for an event-by-event scattering amplitude which contains a noise term. If one starts with a 

I fixed initial condition r(r, Iq), a single realization of the noise leads to an amplitude for a single event, while different 
Oh, realizations of the noise result in a dispersion of the solutions. 

^ It has been observed that, event-by-event, the travelling-wave property is preserved and that the major effect of 

. . , this noise term is to introduce dispersion in the saturation scale Qs from one event to another. The saturation scale 
^ ■ becomes thus a random variable. Very recently, the cumulants of the probability distribution of the saturation scale 
(or, more precisely, of its logarithm), have been computed j^. In the present work, we reconstruct the probability 
distribution from the cumulants and study its different asymptotic behaviours, i.e. the probability for fluctuations 
\ far above, around, or far below the average saturation scale Qg. 

We prove that the probability distribution is Gaussian within a large window around the average saturation scale. 
We also compute the probability distribution outside of this window and compare our analytical predictions with the 
numerical simulations obtained in 0. We justify the use of the Gaussian law previously considered in the literature 
[M 1^, ITol lllf and responsible for a new scaling law of the dipole amplitude: when computing the physical amplitude 
by averaging over all events using a Gaussian distribution, one gets T(r, Y) = T(\og[r'^Ql{Y)\/ ^/Y). 

The structure of the paper is as follows. We start in section by introducing the QCD evolution equations and 
their link with the sF-KPP equation. In section IIIII we compute the probability distribution from the cumulants 
obtained in Q and compare those predictions with numerical simulations. Finally, we discuss the implications on 
physical amplitudes in section Hvl We conclude in section Ivl 
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FIG. 1: These diagrams show the contributions to the evolution of T^^\ the scattering amplitude for a set of two dipoles off 
an arbitrary target represented here by a gray disk. From left to right: a linear BFKL term proportional to T'^', a merging 
term proportional to T*-^' responsible for saturation, and a splitting term proportional to T responsible for fluctuations. 



II. STOCHASTICITY IN HIGH-ENERGY QCD EVOLUTION 



To begin with, let us recall the important points concerning the QCD evolution equations towards high energy. 
To the present knowledge (in the leading logarithmic approximation and in the large- A^c limit), those equations give 
the energy evolution for the scattering amplitude ^T'^'^^) between a system of k dipoles of transverse coordinates 
(xi,yi), . . . , (x/j,yfc) and a generic target. They include three different types of contributions as shown in Fig. 
n BFKL ladders, gluon mergings responsible for saturation corrections, and gluon splitting corresponding to gluon- 
number fluctuations. Summing these contributions results in a infinite hierarchy of equations, containing Pomeron 
loops, in which the evolution of ^T''"')) depends on (T^'^^), (^T^'^^^^) and (T*^'^~^)), coming respectively from the three 
types of contributions discussed previously. 

The solutions of this hierarchy are not yet known and it is useful to deal with a simplified version of it. If we 
integrate out the impact-parameter dependence by performing a coarse-graining approximation 0, we can Fourier 
transform the set of equations to momentum space (r k). The resulting hierarchy becomes equivalent to the 
following Langevin equation 

dyT{k,Y)=a r ^ 
Jo P 



p'^T{p,Y) - k^T{k,Y) k^T{k,Y) 



aT^{k, Y) + a^/2KalT{k, Y) v{k, F), (1) 



where v{k, Y) is a Gaussian white noise satisfying 

{v{k, Y)) = and (i/(fc, Y)v{k' , Y')) = —5 (log(fcVfc'^)) S(Y - Y'). (2) 

an 

The factor k present in front of the noise term comes from a local-noise approximation performed during the coarse 
graining. 

• The linear part of equation is the BFKL equation Its solution Tiin(fc,F) is a superposition of waves: 

TUk,Y)= ^ e-^[^-*"(^)^l ro(7) (3) 

where L — log(fc^//cQ) and Tq{'j) specifies the initial condition, characterised by the reference scale fco- Each 
wave has a different speed v{'y) given by 

vh) = ^ with x(7) = 2V(1) - V(7) - ^(1 - 7) , (4) 

7 

i/;(7) = ^ log[r(7)] being the digamma function. 

• The equation Q without the noise term is the BK equation. It has been shown U to lie in the same universality 
class as the F-KPP equation Q which allowed analytical predictions. The non-linear term damps the growth 
of the amplitude predicted by the linear BFKL equation in such a way that the asymptotic solution Ttw{k, Y) 
of the BK equation is a travelling wave of minimal speed Ttw [L — aVcY) . The minimal speed Vc is obtained for 
a value of 7 that we shall denote 7^ : 

Vc — min v("f) = w(7c) with 7c solution of — x'il) ■ (5) 

7 7 
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For the BFKL kernel one has 7c = 0.6275 and Vc = 4.883. From the position of the travelhng wave one 
obtains the saturation scale Qs{Y), the momentum for which T^^ is constant: 

Ttw(L - av.Y) e--=(^— = ''^ with QUY) = e^^^ . (6) 

• The complete equation is a stochastic equation: it is the BK equation supplemented with a noise term which 
accounts for gluon-number fluctuations. It has been shown that, in the diffusive approximation where the 
integration in is expressed as a second-order differential operator, Eq. lies 's'l in the same universality 
class as the sF-KPP equation 

dtu{x, t) = dlu{x, t) + u{x, t) — u^{x, t) + ^/2Xu{x,t)[^ - u{x,t)]r]{x, t), (7) 

with 

{rj{x, t)) = and (ry(a;, t)Tj{x',t')) = 5{t - t')d{x - x'). (8) 

In the analogy between and aY plays the role of time t while the space variable x is related to L and 
A cx nal . The main analytical predictions concerning ^ are obtained from the study of the sF-KPP equation in 
the weak-noise limit. First, if one compares solutions of ^ with solutions of the BK equation, one obtains that 
each individual event is a travelling wave with a speed smaller than the one predicted from the BK equation. 
Then, if one consider a whole set of events, the main effect of the noise is to introduce dispersion in the position 
of different events. Physically, this means that the saturation scale fluctuates from one event to another and 
we can show, analytically (in the weak or strong noise limit) and numerically, that the average saturation scale 
Os(^) grows as exp(ai;y) with v < Vc and that the dispersion oi\og{Ql{Y)) around that average value increases 
like Vy as expected from a random-walk process. 



III. THE PROBABILITY DISTRIBUTION OF THE STOCHASTIC SATURATION SCALE 

The saturation scale being a stochastic variable, it is characterised by a probability distribution. Let us call P{ps) 
the probability distribution for the variable ps = \og{Ql{Y) j k"^) . From general arguments it was first inferred that 
P{ps) should be Gaussian with an average value {ps) = \og{Ql{Y)/k'^) = avY and a variance — aDY, with v and 
D are coefficients characterising the speed and dispersion of the travelling wavefronts. This was then used extensively 
in the literature JLOj HJ- Recently, the cumulants k„ of the distribution have been computed. The first cumulant 
Ki = (ps) is the average value of ps and the second cumulant K2 = is the variance of the distribution. In the 
Gaussian case higher-order cumulants are zero. The calculation of showed that this was not the case: they find 
that higher-order cumulants arc proportional to the second cumulant. To summarise, one has 

Ki = (ps) = avY , K2 = cr^ = aDY , k„ = ' „ , (9) 

where C(n) is the Riemann Zeta function. Note that these results are valid in the high-energy limit 7cCr^ 1. 
In the weak-noise {a^ ^ 1) limit in which the results were obtained, one has 

.^.(7.)--T^"-^ A„31og[log(l/a-)]A ^^jHfij^^ 

' ''' log2(l/«2) V log(l/«^) ; 31og3(l/a2) ^ ^ 

But the fact that (ps) and are proportional to Y is believed to be more general. This has been shown by an 
analytical study of the strong noise limit and it is confirmed by numerical simulations IT^ for arbitrary values 
of the noise strength (see also fl^ for related numerical studies). One observes that, when the noise strength increases, 
the speed of the wave v decreases and the dispersion coefficient D increases. In what follows, we shall therefore keep 
V and D as parameters. 



A. Analytical results 

In this section, we compute the probability distribution P{ps) from the cumulants. Our starting point is the 
generating function for the moments of P{ps) 

/oo 
dp.e'^P^ P{ps) . (11) 
'OO 
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The cumulants generating function then reads 



C(n)A" 



n>0 



If we use the foUowing integral representation of the Zcta function 



1 



du ■ 



1' 



(12) 



(13) 



r(n) J, 

we are able to compute analytically the cumulants generating function H12|l 

TT^ 7c Jo e" - 1 TT^ 

where 7^; « 0.577216 is the Euler constant. We can then inverse the Laplace transform in (|11|) to obtain the probability 
distribution (c < 7c) 



7_E + i/" 1 

7c 



(14) 



PiPs) = 



c-\-ioo 



dX 

2iTT 



where we have introduced 



exp < — Az — 6A 



q: 



( ^ 



6 = 



37c«T^ 



(15) 



(16) 



z is the distance of the saturation scale to its average value, and 6 is a convenient redefinition of the variance. 
To evaluate further the probability distribution P{ps), we shall perform the integration over A in the saddle point 
approximation. One finds that the saddle point A = 7c(l — e) has to satisfy 



7£;+V(e)-(l-e)^<')(e) 







(17) 



where ip^"''{x) is the polygamma function defined as ■j-^i}j{x). The probability distribution is then given by 



P{Ps 



12^2 



-1/2 



exp 



3cr 



:^7c^(l-e)2^W(£) 



(18) 



Although Eq. (|17l) has no exact analytical solution, one can solve it in three interesting limits: {i) z/b ^ which 
corresponds to e ^ 1; (ii) z/b — > 00 related to the limit £ — > 0; and (Hi) z/b — > —00 equivalent to the limit e — *■ 00. 
We detail these limits hereafter. 



(i) When z/b ^ or equivalently \z\ <^ "ic'^'^ ■ Eq. H17|l becomes 



7i=; + V(l)-2(l-e>W(l) 



= 



e = 1 - 



iz 



where we have used ■0(1) = —1e and 0^^'(1) = t:"^ After replacement in H18() we get 



P{Ps) 



1 



V2^ 



1 



9C(3) z 



exp 



z 

'2^ 



We kept the first subleading term to show that this distribution has a maximum for the constant value 

. -9C(3) 



-1.75. 



(19) 



(20) 



(21) 



The most probable value for the saturation scale is therefore not the average value {ps)- For asymptotically 
large energies, we recover 



P{Ps) 



: exp 



z 

2^ 



(22) 



The probability distribution of fronts is Gaussian around the average. We shall discuss later on the fact that 
this is the dominant behaviour at high-energy which justifies the use of a Gaussian distribution in |5llR ll0llll[ . 
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(ii) When z/b ^ +00 or equivalently z ^ jcO''^. We are sensitive to small values of e. One can thus simplify 
Eq. l(T7|) and obtain 



TT 







3 7ctT2 



1 



Teg 
4z 



After a bit of algebra, this turns into the following limit for the probability distribution 

,2 \ 3/4 



3I/4 



exp 



3 7cf72 



This corresponds to a power-law tail P{ps) ^ (Qs/Qs) ^''^ at very large Qs- 
(iii) When z/b — > —00 or equivalently z ^ — 7cCr^. We have e ^ 00, and Eq. H17|l turns into 

1 



T+1E + l0g(f ) + 1 . 

e 



The probability (|18|1 is then found to be 







exp 



PiPs) 



TT 1 



6 cr 



exp 



1 + 7B 372a2 



exp 



3 7cCr2 



3 700"^ 



1 



7b 



1. 



7e 



(23) 



(24) 



(25) 



(26) 



6 7cCT^ 2 TT' 

This is a Gumbel distribution, which goes to zero very fast when z ^ — 7cCr^. 

As a function of z starting from z — —00, the transition between the regime H26|l and the regime H22(l happens for 
z = — 7cO'^ and the transition between the regime H22|) and the regime (|24|l happens for z = 7cCr^. For those two points, 
Ps — (ps) i 7cC^, the probability is of order e^'^'^'^ ja which is very small. This means that the probability P{ps) 
is not Gaussian only for very improbable fluctuations. To describe the bulk of the events, the Gaussian distribution 
is a good approximation. 



B. Numerical results 



In this section, we shall compare the analytical prediction H18|l derived in the previous section with the numerical 
simulation of introduced in 9J . To fix things properly we need first to recall a few points from . In that study, 
we start with a fixed initial condition and evolve numerically the QCD Langevin equation with different realisation 
of the noise term (the precise method used to solve the equation is not important for our purposes here, so we refer 
to 01 for details). The momentum k is discretised in bins of L = log(fc^/fcg) between Lmin and Lmax- Hence, the 
numerical simulations of results in a set of events Ti(Lj, Y), with 2 = 1,..., ricv, for each event i, we get the 
rapidity evolution of the amplitude in each momentum bin. 

As expected from the analytical study of |Q in the weak-noise limit, the numerical studies have confirmed that 
each event is a travelling wave whose speed decreases when the strength of the fluctuations increases. Also, at large 
rapidities, the position of the wavefront shows a dispersion proportional to \/Y if one consider a whole set of events. 
Practically, those results can be summarised by saying that, for the ith event, we have 

T,{L,Y) = T,{L-L,^,(Y)) with = = log(Q2(y)/fc2)_ (27) 

The position Lsj{Y) of the wavefront for one event can be extracted from the numerical simulations by solving 
Tj(i, Y) = Tq for a fixed Tq (here, we adopted Tq — 0.2) and for different values of the rapidity. By measuring (the 
logarithm of) the saturation scale, one can obtain its statistical distribution. More precisely, for each rapidity, we 
construct a histogram by counting, for a set of events, the number of event for which the saturation scale is in each 
momentum bin. 

In order to compare with the analytical predictions, we shall, instead of the histogram for ps, use the distribution 
for z = ps — (Ps)- Normalised in such a way that its integral is 1, it is directly comparable with H18|l if one fixes the 
dispersion parameter. For the latter we adopt the value computed numerically from the histogram of ps or z. In 
other words, we fix the two lowest-order cumulants (ps) and a^, and compare the resulting analytical and numerical 
probability distribution at different rapidities. 
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FIG. 2: Comparison of the probability distribution 1)18^ as a function of z = — (pa) with resuhs from numerical simulations 
with a = 1 (top) and k = 5 (bottom). The comparisons are made for Y — 20, 30, 40, from left to right. 

The results are shown in Fig. 13 for two sets of 10000 events taken from 0. Those two sets corresponds to 
different values of the noise strength k = 1Qttk/N^ = 1 and 5, the value of a being fixed to 0.2. For both cases, 
we show the probability distribution for rapidities Y = 20, 30 and 40. Those rapidities are in the region where the 
travelling wavefront has reached its asymptotic behaviour and the dispersion of the wavefronts is proportional to Vy 
as expected. 

It is obvious from Fig. |21that the agreement is excellent with predictions from IjlSI) . One can see that, in the bulk 
of the distribution, the probability is Gaussian, up to a shift of the maximum to a negative value which is consistent 
with (|21|l . We also observe that the probability falls very fast in the infrared and has a tail which favours fluctuations 
to large values of the saturation momentum. However, the deviations from a Gaussian behaviour only appear for 
events with a very small probability. 



IV. CONSEQUENCES FOR THE PHYSICAL SCATTERING AMPLITUDE 



We shall now compute the behaviour of the average scattering amplitudes {T{p,Y)) where p = log[l/(fcQr^)], using 
the probability distribution obtained in the previous section. For simplicity, we consider the amplitude in coordinate 
space. The average amplitude can be directly computed from the probability P(ps) using 



dPsP{Ps)T{p- ps), 



(28) 



where T{p—ps) describe the event-by-event amplitude which features the travelling-wave property w.r.t. the saturation 
scale Ps- Inserting the probability distribution (jlSj) into H28|l and using the McUin representation 



T{p) 



dX 
2i7r 



T(A), 



we obtain for the physical amplitude 



{T{p,Y)) = J ^nX) exp|-AZ-6A 



7b + V' 1 

7c 



(29) 



(30) 
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We have introduced the convenient variable 

Z = p-{ps) = -log {v^QliY)) 



(31) 



Since the function T[\) is slowly varying the remaining integration will basically be sensitive to the same saddle 
point H17|l as for the computation of the probability. The main difference lies in the fact that, because of the unitarity 
constraint, r(A) has a pole in 7 = 0. Hence, in order to get relevant expressions in the different physical limits, the 
following assumption for the event-by-event wavefront is sufficient: 



1 P< Ps 

exp[-7o(p-Ps)] P>Ps 



f(A) 



1 
A 



70 



(32) 



This satisfies both the travelling-wave property ® and unitarity requirements. The exponent 70 can be chosen-'^ 
between 7^ and 1. 

We can now introduce (|32l) in H3U|) and we get 



{T{P,Y)) = 



7o 



dX 



2«7r A(7o — A) 



exp < —\Z — bX 



7_B -I- ?/; 1 - 



7c 



(33) 



where the contour of integration is now restricted to values of c such that < c < jc- The remaining integration over 
A can be computed in the saddle-point approximation and the condition H17|) is still appropriate for H33|) . Again, the 
physical picture arises when we consider the three different limits |.^| ^ Ic'^'^, Z 3> 7cO'^ and Z ^ — 7ccr^, which we 
detail hereafter. 

(z) Case Z 3> 7cC^: as for the computation of the probability, the saddle point is obtained by expansion around 
e ^ or A — !■ 7c and the same result is obtained, up to an extra factor (70 — A)^^. For 70 > 7c, this leads to 



{np,Y)) 



31/4 



27rcr(7o - 7c) V Z 
If 7o — 7ci we get an additional contribution from the pole^ 



3/4 



cxp(-7cZ). 



{np,Y)) 



27cCr 



7cCT" 



1/4 



exp(-7cZ). 



(34) 



(35) 



{ii) Case Z <C — 7cCr^: this case is slightly more intricate. Indeed, since A —00, the contour of integration must 
first be moved towards c = A. The integration therefore receives a contribution of 1 coming from the pole at 
A = 0. The final result reads 



{T{p,Y))^l 



TT 1 

6 7cCT 



exp 



Z 



2 7cCT^ 



+ 7:(l+7£;) 



37cV2 



exp 



3 7ccr2 



1 



7s 



(36) 



[iii) Case \Z\ <C 7cO'^: for that last case, we have to distinguish between two possibilities: Z > and Z < 0. Let 
us start with the positive values of Z: as for the computation of the probability, this situation is sensitive to 
A — > 0. Hence, we need to properly take into account the prefactor A^^ = exp[— log(A)]. This gives an additional 
contribution — A~^ to the saddle point equation which becomes 



Z-cr^A 



1 

= 

A 



A = 



2a2 



4^2 



(37) 



^ Although we keep it as a variable, one should probably adopt 70 = 7c. This exponent in the tail differs from the usual 70 = 1, however, 
since the results are obtained in the weak-noise limit, the tail of the wavefront extends far above the saturation momentum and the 
proper matching with the exp{— p) behaviour should only lead to higher-order corrections. 

^ If we treat properly the pole inside the saddle-point equation, it only introduces a factor (7c — A)~^ which is subleading compared to 
the (7c — A)~^ from the ^^^^ function. Hence the saddle point remains unchanged. 
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FIG. 3: Comparison of the physical amplitude (T{p,Y)) plotted as a function of Z/a and the error function obtained as the 
asymptotic behaviour. 

If, in addition, one requires Z ^ a, we recover X = Z/a^ as for the computation of the probability. The average 
amplitude is then the same as the probability, up to an extra factor A^^, which gives 



For the case where Z is negative, we first move the integration contour on the left side of the pole at A = 
which gives a contribution of 1 and the remaining integration is computed in the same way as for positive Z. 
That is, for — 7ca'^ Z —a, we obtain 



where we have explicitly emphasised the sign of the second term. The two results ()38|) and H39fl can be sum- 
marised in only one formula: 



where erfc(a;) is the complementary error function. 

The results (|35ll . (|3t)|) and (|40|l call for a number of comments. The most important one is that (for very high 
energies such that 7cO'^ 3> 1) one recovers the error function as expected previously in the literature 0, 0, liol 
and the scattering amplitude satisfies diffusive scaling as it depends only on the ratio Z / <t. This result is valid in a 
very large window around the average saturation scale {Z ^ TcC^). Since this can directly be derived from the front 
(|32|l convoluted with the Gaussian probability distribution H22|l . this means that it is sufficient to consider a Gaussian 
distribution in order to compute the scattering amplitudes at high-energy. In addition, the error function fully comes 
from yO < ps in H32() , which proves the expected result [Tol ITU that the amplitude is dominated, up to very large 
values of p, by fronts at saturation i.e. by hlack spots. Finally, let us notice that our argument does not depend on 
the particular choice for the event-by-event front. Indeed, the pole in A~^ comes from the saturated part of (|32|l and 
thus the condition that the event-by-event amplitude satisfies unitarity is sufficient to obtain (|40|) . 

To illustrate that behaviour, we have computed from H28() the amplitude obtained using an event-by-event front 
0(ps — p) {Q{x) being the Heaviside function). We plot on Fig. Othe result for different values of ct as a function 
of the high-energy scaling variable Z/a. We clearly observe that when a increases {i.e. when energy increases), we 
converge to the universal error function (|40|l . 

The behaviour in the far tail Z ^ 7cO''^ is also quite interesting. Indeed, if we adopt 70 = 7c, the amplitude 
(|35|l decreases as exp(— 7c2') as expected, but it receives contributions both from the event-by-event amplitude and 
from the probability distribution. Finally, deep inside the saturation region we obtain that the amplitude reaches the 
unitarity limit like 1 — ci exp[— C2(r(5s)^'^^] where ci, C2 and C3 can be read on Eq. (|36|) . 




(38) 




(39) 




(40) 



9 



V. CONCLUSION AND DISCUSSION 



Let us now summarise the results obtained in this paper. First, we have shown that it is possible to compute the 
probability distribution for the position of the wavefront satisfying the stochastic F-KPP equation, which corresponds 
to the saturation scale in QCD. Our starting point is the cumulants of this distribution computed previously in Q. The 
probability distribution is computed in a saddle point approximation. This allows us to compute its leading behaviour 
in the vicinity of the average saturation scale (p^) where the probability is Gaussian, as well as for improbable events 
far above or below (ps)- 

Then, we have checked that the predictions for the probability distribution are in good agreement with the results 
obtained from the numerical simulations of the QCD Langcvin equation , once one fixes the average position and 
its dispersion from their numerical values. Since the numerical simulations are not obtained in the weak noise limit 
and not directly for the sF-KPP equation, this proves once again that the probability distribution derived in this 
paper has some universal properties. 

From the probability distribution we are able to deduce the physical scattering amplitude. Again, we have considered 
the same interesting asymptotic behaviours as for the probability distribution: the vicinity of the average saturation 
scale, deep inside the saturation regime or in the dilute domain. Within this analysis, the most important result is that 
at high-energy the probability distribution can be considered as Gaussian over an very large domain |z| <C 7cC^ cx Y. 
As a consequence the scattering amplitudes are dominated by events which are at saturation. 

This validates the approach adopted in Refs. 0,0], and leads to the fact that scattering amplitudes scale as a 
function of log[r^Qs(^)]/v^: ^ property which is known as diffusive scaling and may have important consequences 
on LHC physics. The domain over which the probability distribution may be considered as Gaussian extends over 
a region of ps which satisfies \ps — {ps)\ ^ TcC^. As the dispersion a grows like VY, this domain becomes larger 
and larger with increasing rpidity. The same argument holds for the region \Z\ ^ 7cO'^ for which the scattering 
amplitude satisfies diffusive scaling. Those facts also corroborate the result from J,3] that, in the strong-fluctuation 
limit, the scattering amplitude is given by the error function H4U|) . which is a superposition of Heaviside functions 
with a Gaussian distribution. 

Of course, the analytical behaviour of the cumulants associated with the distribution of Ps for any value of a^. and 
in the case of the full QCD equation Q still remains an open question. However, our analytical and numerical results 
show that they are compatible with a universal probability distribution. This deserves more detailed studies. 
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